#include <iostream>
#include <fstream>
#include <vector>
#include <string>

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#include<iostream>
#include<fstream>
#include<typeinfo>

#include "util.h"
#include "diffuse.h"

// 定义方阵的维度
#define nMAX 10  

using namespace std;


int main()
{

	char str1[64], str2[64], str3[64];

	sprintf(str1, "output_data");
	if (!access(str1, 0) == 0) mkdir(str1, 0777); // if  no "data" folder, construct it for output date

	char   buf_data[40];

	int ndim = 3;
 	
	if (ndim > nMAX) { cout << "矩阵阶数超过最大数nMAX !!!" << endl;  exit(-1); }

	// double  a[nMAX][nMAX] = { 0 };//定义一个nMAX*nMAX*的矩阵，用于存放数据
	double  **a = dmatrix(1, ndim, 1, ndim); // 或者我们可以动态开辟数组，这种方式记得在程序结束后释放内存。使用方式 free_dmatrix(a, 0, nMAX, 0, nMAX);
	
	// 怎样读取数据文件，比如data.txt中存储的是矩阵A数值
	ifstream infile;//定义读取文件流，相对于程序来说是in
	infile.open("data.txt");//打开文件
	for (int i = 1; i <= ndim; i++)//定义行循环
	{
		for (int j = 0; j <= ndim; j++)//定义列循环
		{
			infile >> a[i][j];//读取一个值（空格、制表符、换行隔开）就写入到矩阵中，行列不断循环进行
		}
	}
	infile.close();//读取完成之后关闭文件


	ifstream infile2;
	infile2.open("data_b.txt");		// 打开文件
	for (int i = 1; i <= nMAX; i++) // 定义行循环
	{
		infile2 >> b[i]; // 读取一个值（空格、制表符、换行隔开）就写入到矩阵中，行列不断循环进行
	}
	infile2.close(); // 读取完成之后关闭文件


////////////////////////////////////////////////////////////

	// Cholesky分解
	double  **a = dmatrix(1, ndim, 1, ndim); // 开辟内存，用存储LU分解后的下三角形矩阵L
	double  **L = dmatrix(1, ndim, 1, ndim); // 开辟内存，用存储LU分解后的上三角形矩阵U

	/*
	对矩阵a赋值，
	Cholesky分解后的矩阵存储在L矩阵中
	*/
	sprintf(buf_data, "output_data");
	


  // 动态开辟的内存，程序结束前，最好释放掉， 不做内存释放，程序运行后，系统也会自动释放内存
   free_dmatrix(a, 1, ndim, 1, ndim);
   free_dmatrix(L, 1, ndim, 1, ndim);
 
   return 0;
}



void Thomas(doubl **a, double *x, double *f, int n)
{
    // ax=f , 追赶法求解三对角线性方程组, 算法2.11
	// a[3][1,..n]
	// x[1,..n]
	// f[1,..n]
	// n: number of unknown

	//a[1]: 三对角中的平行对角线下面的一行, 算法中的a_i
	//a[2]: 三对角中对角线行, 算法中的b_i
	//a[3]: 三对角中的平行对角线上面的一行, 算法中的c_i

	double *apha = dvector(1, n);
	double *beta = dvector(1, n);
	double *y    = dvector(1, n);

	apha[1] = a[2][1];          
	beta[1] = a[3][1]/a[2][1];  // b  
	y[1]    = f[1]/a[2][1];

	for (int i = 2; i <= n - 1; i++)
	{
		apha[i] = a[2][i] - apha[i - 1] * beta[i - 1];
		beta[i] = a[3][i] / apha[i];
		y[i] = (f[i] - a[1][i - 1] * y[i - 1]) / apha[i];
	}

	apha[n] = a[2][i] - apha[n - 1] * beta[n - 1];
	y[n] = (f[n] - a[1][n - 1] * y[n - 1]) / apha[n];
	x[n] = y[n];

	for (int i = n-1; i >= 1; i--)
	{
		x[i] = y[i] - beta[i] * x[i + 1];
	}

	// 计算函数中动态申请的内存，在算法计算使用完成后，释放掉，防止多次调用后，内存泄漏
	free_dvector(apha, 1, n);
	free_dvector(beta, 1, n);
	free_dvector(y,    1, n);

	return;
}


void print_vector(double *x, int n, char* buf)
{
	char buffer_x[32];
	
	FILE *fx;

	sprintf(buffer_x, "./%s/x.m", buf);

	fx = fopen(buffer_x, "w");
	for (int i = 1; i <= n; i++)
	{
		fprintf(fx, " %16.14f ", x[i]);
		fprintf(fx, "\n");
	}
	fclose(fx);

	return;
}



void print_matrix(double** phi, int n, char* buf)
{

	char buffer_phi[32];

	FILE *fphi;

	sprintf(buffer_phi, "./%s/a.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	return;
}


void print_matrix2(double** phi, double** eta, int n, char* buf)
{

	char buffer_phi[32], buffer_eta[32];

	FILE *fphi;
	FILE *feta;


	sprintf(buffer_phi, "./%s/a.m", buf);
	sprintf(buffer_eta, "./%s/L.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	feta = fopen(buffer_eta, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(feta, " %16.14f ", eta[i][j]);
		}
		fprintf(feta, "\n");
	}
	fclose(feta);

	return;
}



void print_matrix3(double** phi, double** eta, double **xx,  int n,  char* buf)
{

	char buffer_phi[32], buffer_eta[32], buffer_xx[32];

	FILE *fphi;
	FILE *feta;
	FILE *fxx;

	sprintf(buffer_phi, "./%s/a.m", buf);
	sprintf(buffer_eta, "./%s/L.m", buf);
	sprintf(buffer_xx,  "./%s/U.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	feta = fopen(buffer_eta, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(feta, " %16.14f ", eta[i][j]);
		}
		fprintf(feta, "\n");
	}
	fclose(feta);


	fxx = fopen(buffer_xx, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fxx, " %16.14f ", xx[i][j]);
		}
		fprintf(fxx, "\n");
	}
	fclose(fxx);

	return;
}


